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ABSTRACT Strong evidence supports the idea that fatty acids rather than carbohydrates are the main energy source of Mycobacte- 
rium tuberculosis during infection and latency. Despite that important role, a complete scenario of the bacterium's metaboUsm 
when lipids are the main energy source is still lacking. Here we report the development of an in vitro model to analyze adapta- 
tion of M. tuberculosis during assimilation of long-chain fatty acids as sole carbon sources. The global Upid transcriptome re- 
vealed a shift toward the glyoxylate cycle, the overexpression of main regulators whiB3, dosR, and RvOOSl, and the increased ex- 
pression of several genes related to reductive stress. Our evidence showed that Upid storage seems to be the selected mechanism 
used by M. tuberculosis to ameliorate the assumed damage of reductive stress and that concomitantly the bacilli acquired a 
slowed-growth and drug-tolerant phenotype, all characteristics previously associated with the dormant stage. Additionally, in- 
tergenic regions were also detected, including the unexpected upregulation of tRNAs that suggest a new role for these molecules 
in the acquisition of a drug- tolerant phenotype by dormant bacilli. Finally, a set of lipid signature genes for the adaptation pro- 
cess was also identified. This in vitro model represents a suitable condition to illustrate the participation of reductive stress in 
drugs' activity against dormant bacilli, an aspect scarcely investigated to date. This approach provides a new perspective to the 
understanding of latent infection and suggests the participation of previously undetected molecules. 

IMPORTANCE Mycobacterium tuberculosis establishes long-lasting highly prevalent infection inside the human body, called latent 
tuberculosis. The known involvement of fatty acids is changing our understanding of that silent infection; however, question of 
how tubercle bacilli globally adapt to a lipid-enriched environment is still an unanswered. With the single change of providing 
fatty acids as carbon sources, the bacilli switch on their program related to dormant stage: slowed growth, accumulation of lipid 
bodies, and development of drug tolerance. In this stage, unexpected and previously unknown participants were found to play 
putatively important roles during the process. For the first time, this work compares the global transcriptomics of bacteria by 
using strand-specific RNA sequencing under two different growth conditions. This study suggests novel targets for the control of 
tuberculosis and provides a new straightforward in vitro model that could help to test the activity of drugs against dormant ba- 
ciUi from a novel perspective. 



Received 4 April 2014 Accepted 8 April 2014 Published 20 May 2014 

Citation Rodriguez JG, Hernandez AC, Helguera-Repetto C, Aguilar Ayala D, Guadarrama-Medina R, Anzola JIVl, Bustos JR, Zambrano IVIM, Gonzalez-y-Merchand J, Garcia MJ, 
Del Portillo P. 2014. Global adaptation to a lipid environment triggers the dormancy-related phenotype of Mycobacterium tuberculosis. mBio 5{3):e01 125-14. doiilO.l 128/ 
mBio.Ol 125-14. 

Editor Carol Nacy, Sequella, Inc. 

Copyrigiit © 2014 Rodriguez etal. This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-ShareAlike 3.0 Unported 
license, vjhlch permits unrestricted noncommercial use, distribution, and reproduction in any medium, provided the original author and source are credited. 
Address correspondence to Patricia Del Portillo, pdelportillo@corpogen.org, or Maria J. Garcia, marlaj.garcia@uam.es. 



It is estimated that 2 billion people are latently infected with 
Mycobacterium tuberculosis, representing potential cases of reac- 
tivation and transmission and hampering better control of the 
disease (1). As tuberculosis (TB) infection progresses, the bacilli 
face adverse environmental conditions to which they respond by 
entering into a silent phase known as dormancy. In this phase, the 
bacteria are able to survive for the entire life of the host causing 
latent infection. Eventually, they can also reactivate and cause ac- 
tive disease when the host's immune system becomes weakened 
(2). Bacterial adaptation to the different stages of the infection 
involves fine regulation of gene expression. Global transcriptome 
analysis has identified critical regulators of gene expression such as 
the p?ioP virulence regulator (3), the kstR cholesterol regulator (4), 



the dosR hypoxia regulator (5), or RvOOSl, which was also recently 
proposed as a hypoxia regulator (6). The identification of noncod- 
ing RNAs in M. tuberculosis also suggests that these molecules 
could regulate microbial adaptation (7, 8). During latency, M. tu- 
berculosis is presumed to reside in lipid-rich foamy macrophages 
(9) and accumulate internal lipid droplets full of triacylglycerols 
(TAGs) that are used subsequently as an energy source for persis- 
tence (10). Nevertheless, the bacilli continue replication (11) and 
retain low but detectable metabolic activity (12). It has also been 
shown that M. tuberculosis uses host lipids, in particular, fatty 
acids and cholesterol, as energy sources during intracellular 
growth and persistence (13-15). 

Although several in vitro models have been used to mimic the 
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TABLE 1 Transcriptomic profile of M. tuberculosis in different carbon sources'^ 




No. (%) of reads in following growth phase: 






Reads 


DE 


DS 




ro 


All 


13.12 


19.29 


18.314 


16.054 


All mapped 


8.17 (62.3) 


10.78 (55.9) 


9.56 (53.4) 


8.01 (51.2) 


Mapped without rRNA 










Mapped to CDS 


0.52 (71.3) 


0.51 (57.6) 


0.28 (61.6) 


0.10 (13.4) 


Mapped to ICRs 


0.17 (22.4) 


0.25 (28.0) 


0.149 (32.9) 


0.62 (83.0) 


Mapped to antisense CDS 


0.05 (6.4) 


0.13 (14.4) 


0.02 (5.4) 


0.03 (3.6) 



" Total numbers (in millions) of reads determined under the different conditions tested (upper part) and numbers (percentages) of reads determined by excluding those 
corresponding to rRNA (lower part) are shown. 



environment of dormant bacilli during latency (16-18), cultures 
grown in the presence of lipids, the most abundant molecules 
surrounding bacilli in vivo, have not been extensively analyzed. 

To gain insight into the process of M. tuberculosis adaptation to 
a fatty acid environment, we developed an in vitro model where 
bacilli were grown in even-length long-chain fatty acids (LC-FAs) 
as the sole carbon source. Even-length LC-FAs were selected be- 
cause they are the main LC-FAs found in human cells (19) and are 
a major component of TAGs (20). Although other in vivo lipid 
energy sources cannot be ignored, it is possible that for long-term 
infection, bacilli might benefit by selecting even-length LC-FAs, 
avoiding the toxicity associated with other lipid sources (21). 

We analyzed changes in the global transcriptome of M. tuber- 
culosis by using high-throughput methods such as strand-specific 
RNA sequencing (ss-RNA-seq) (22). 

RESULTS AND DISCUSSION 

The global transcriptome of M. tuberculosis in the LC-FA in 
vitro model. In order to mimic the nutrient conditions presumed 
to be encountered by M. tuberculosis during in vivo infection (23), 
we cultured M. tuberculosis H37Rv in a medium supplemented 
with a mixture of even-length LC-FAs as the sole carbon source. 
Growth in LC-FAs was compared to growth in dextrose as a con- 
trol. We found that growth in LC-FAs was slower than growth in 
dextrose but reached similar saturation points (see Fig. SI in the 
supplemental material). 

RNA was analyzed by ss-RNA-seq with material from bacilli 
grown in dextrose to the exponential (DE) and stationary (DS) 
phases and in LC-FA to the exponential (FE) and stationary (FS) 
phases (see Fig. SI). The sequence depth obtained showed good 
coverage of the M. tuberculosis genome, as indicated by the pla- 
teaus reached in all instances by mapping reads against coding 
sequences (CDS) (see Fig. S2A). A normal distribution of the data 
was achieved in all cases (see Fig. S2B) (24). The number of se- 
quences retained after processing for quality was 13.1 to 19.3 mil- 
lion reads (Table 1). To facilitate the comparison of gene expres- 
sion, the data were normalized as reads per kilobase per million 
reads (RPKM). Only genes showing >1 RPKM were considered 
for analysis. Statistical analysis for significant differential gene ex- 
pression was performed with Fisher exact test (see Table SI A and 
B in the supplemental material). That analysis gave results consis- 
tent with those obtained with the DEGseq package for biological 
replicates (see Table SIC and D) (see Materials and Methods). 

We identified all 45 tRNAs (25) (see Table S2A), all of the most 
highly expressed small RNAs (7), and most of the previously iden- 
tified noncoding RNAs (8) (see Table S2B). During growth in 



LC-FAs, gene expression (reads mapped to CDS) was greatly re- 
duced compared with the control, particularly in the stationary 
phase (Table 1). This low gene expression contrasts with the high 
number of reads that mapped to intergenic regions (IGRs), sug- 
gesting close control mediated by small RNAs in cultures grown 
with LC-FAs. 

Consistent with the high number of reads mapping to IGRs in 
the FS culture, the four most highly expressed IGRs corresponded 
to this sample (see Table S2B and C). Two of them, namely, 
MTS2823 and MTS0977 (see Table S2C), have also shown in- 
creased expression during infection (7), suggesting that they 
might be important regulators and highlighting the association of 
the conditions used in our model with the conditions encountered 
by the bacilli in vivo. 

More remarkable was the overexpression of tRNAs in FS (see 
Table S2A), which contrasts with the few transcripts mapped to 
CDS in FS and which, presumably, was a reflection of low trans- 
lation activity (Table 1 ) . Aminoacyl tRNAs are involved in cellular 
processes other than protein synthesis, such as the modification of 
bacterial membrane lipids via the multiple-peptide resistance fac- 
tor (MprF) protein (26). MprF attaches positively charged amino 
acids to the polar heads of phosphatidylglycerol, effectively neu- 
tralizing the negative charge of the membrane (27) and thus de- 
creasing membrane permeability and susceptibility to antibacte- 
rial peptides and cationic antibiotics (28). The tRNAs involved in 
these changes are tRNA-Lys, tRNA-Ala, and less frequently tRNA- 
Arg (26, 29), the same tRNAs that were expressed more in FS than 
in DS in our model (Fig. 1). We hypothesize that, in the presence 
of FA, bacilli modify their envelope by aminoacyl- tRNAs, confer- 
ring adaptation to stress and possibly increasing tolerance to an- 
tibiotics, a process associated with the dormant stage of M. tuber- 
culosis (17). Whether this attractive mechanism to modify the cell 
envelope is, in fact, involved in the adaptation and survival of the 
bacilli during long-term infection remains to be investigated. 

Comparative analysis of CDS between growth in dextrose and 
LC-FA showed that majority of genes had similar expression lev- 
els, both at exponential and stationary phases (see Fig. S3A and B 
in the supplemental material), with only 10% of the potential CDS 
being differentially expressed (Fig. 2; see Table S3). This analysis 
identified genes with differential expression in all of the functional 
categories of TubercuList (http://www.tuberculist.epfl.ch), with 
the exception of genes belonging to the PE/PPE category in the 
exponential phase (Fig. 2; see Table S3). On the contrary, the 
PE/PPE category was overrepresented in the stationary phase 
(Fig. 2 and Table S3), suggesting that when M. tuberculosis grows 
in a fatty acid environment, these abundant proteins (25) might 
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FIG 1 Differential expression of tRNAs in tlie fatty acid model. Fold change comparison of tRNAs between the stationary and exponential phases under both 
culture conditions. Fold changes were calculated by calculating DS/DE (black bars) and FS/FE (light brown bars) RPKM ratios. 



play an additional role in antigenic variation and immune evasion 
(30, 31). Contrary to our expectations, genes belonging to the 
functional category "lipid metabolism" were not overrepresented 
in either the exponential or the stationary phase (Fig. 2; see Ta- 
ble S3), a result that could be attributed to the gene redundancy 
present within this category (25). Among the most highly differ- 
entially expressed genes in FS, we identified pcfcA (P = 4"''') and 
tgs-1 {P = 2~'^) (see Table SIF), both ofwhichplaycritical roles in 
the lipid metabolic pathways of M. tuberculosis (10, 32). As hap- 
pens with pckA, the gene icl-1 was also upregulated within the 
trichloroacetic acid (TCA) cycle (see Table SIE and F), consistent 
with the central role of these enzymes in M. tuberculosis lipid me- 
tabolism (13). Because of the use of only even-length LC-FAs, we 
did not expect production of propionyl coenzyme A (CoA) from 
lipid degradation. In agreement with this, we did not detect in- 
creased expression of genes belonging to the methylcitrate and 
methylmalonate pathways (33). 

The functional category of "regulatory proteins" was under- 
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FIG 2 Genes differentially expressed in the fatty acid model. Numbers of genes differentially expressed 
by functional categories according to TubercuList (http://www.tuberculist.epfl.ch/) are shown. Circle 
sizes correspond to the numbers of genes with significant differential expression in fatty acid medium 
versus dextrose medium. (A) Exponential phase. (B) Stationary phase. 



represented in the stationary phase, highlighting the involvement 
of four genes upregulated in the same phase (FS) (Fig. 2; see Ta- 
ble S3). Two of these correspond to the known regulators WhiB3, 
which is associated with redox balance (34), and DosR, which is 
related to redox balance, and a major regulator of hypoxia in 
M. tuberculosis (5). As a result, a large number of DosR- regulated 
genes were also expressed more in FS than in DS (Table 2; see 
Table SIF). We also identified the recently reported regulator 
RvOOSl, which is considered a global hypoxic response regulator 
(6) and the presumptive transcriptional regulator nrdR (http: 
//www.tuberculist.epfl.ch) . 

Given that WhiB3 and DosR, two heme sensor proteins related 
to intracellular redox balance (34, 35), were identified as key com- 
ponents of the adaptation of M. tuberculosis to growth in LC-FA, 
we searched for significantly overexpressed genes involved in re- 
dox balance. Independently of the phase of growth, bacilli cul- 
tured in LC-FAs had a higher percentage of significantly overex- 
pressed genes associated with reductive stress (see Table S4A). 

These results suggest that growing M. tu- 
berculosis in a fatty acid medium puta- 
tively causes a burst in reductive stress 
and thus the required molecules to main- 
tain redox balance would be induced in 
order to allow cell growth under these 
conditions. An imbalance of the redox 
equilibrium leads to either reductive or 
oxidative stress responses, the latter well 
documented to be involved in M. tuber- 
culosis infection (36). However, the con- 
tribution of its counterpart, namely, re- 
ductive stress, has barely been considered 
thus far and occurs when reducing equiv- 
alents, such as NADH, FADH, or NA- 
DPH, are abnormally increased (34, 37). 
In our LC-FA cultures, the reducing 
equivalents that could be generated by 
both jS-oxidation of LC-FAs and TCA 
metabolic activities would result in re- 
ductive stress that must be cleared for 
survival via reductive sinks, such as poly- 
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RPKM ratio'' 




P value' 












dene tag 


name 


categoi'y'^ 


runction (reierence) 


FE/DE 


FS/DS 


CP ire FITh 

rr, vs ur, 


to vs Lfj 


Rv0079 




10 


Unknown; probable transduction regulatory 
protein (42) 


1.59 


2.10 


1.09 E-12 


3.37 E-07 


RvOOSO 




10 




d 


9.00 




1.19 E-09 


RvOOSl 




9 


Global regulator in hypoxia 





1.10 


— 


0.014 


Rv0243 


jadA2 


1 


Lipid degradation 





1.60 


— 


9.87 E-05 


Rv0569 




10 







3.07 


— 


0.90 E-14 


Rv0570 


nrdZ 


2 


Ribonucleoside-diphosphate reductase 





2.67 


— 


0.004 


Rvl436 


gap 


7 




1.15 


— 


0.028 


— 


Rvl652 


argC 


7 


Arginine biosynthesis 


— 


3.00 


— 


8.33 E-05 


Rvl733c 




3 


Unknown; probable conserved transmembrane 
protein 


2.77 


1.50 


2.62 E-08 


0.012 


Rvl738 




10 






1.15 


— 


0.90 E-14 


Rvl813c 




10 




3.53 




1.09 E-12 




Rvl996 




0 


Unknown; universal stress family protein 


1.93 


1.19 


1.09 E-12 


6.27 E-06 


Rv2005c 




0 


Unknown; universal stress family protein 


— 


6.50 


— 


3.15 E-13 


Rv2007c 


fdxA 


7 


Ferredoxin FdxA 





2,18 


— 


1.09 E-12 


Rv2030c 




10 




— 


5.25 


— 


3.88 E-10 


Rv2031c 


hspX 


0 


Stress protein induced in hypoxia 


— 


6.05 


— 


0.90 E-14 


Rv2623 


TB31.7 


0 


Unknown; universal stress family protein 


— 


3.08 


— 


0.90 E-14 


Rv2625c 




3 


Unknown; probable conserved transmembrane 
protein 


— 


3.20 


— 


1.44 E-06 


Rv2626c 


hrp-1 


10 




— 


4.67 


— 


0.90 E-14 


Rv2627c 




10 






2.43 


— 


2.96 E-06 


Rv2628 




10 






1.40 




1.18 E-05 


Rv2629 




10 






1.13 




4.16 E-08 


Rv3127 




10 






1.33 




0.003 


Rv3130c 


tgs-1 


1 


TAG 




2.31 




2.07 E-13 


Rv3131 




10 


Unknown; putative NAD(P)H nitroreductase 


2.89 


4.00 


1.09 E-12 


9.04 E-13 


Rv3133c 


devR 


9 


Regulator part of the Dos TCS" 




2.60 




6.54 E-05 


Rv3134c 




0 


Unknown; universal stress protein family 


1.03 


3.43 


5.99 E-06 


3.19 E-10 



" Functional categories based on http://www.tuberculist.epfl.ch: 0, virulence, detoxification, and adaptation; 1, lipid metabohsm; 2, information pathways; 3, cell wall and cell pro- 
cesses; 7, intermediary metabolism and respiration; 9, regulatory proteins; 10, conserved hypothetical proteins. 

^ RPKM ratio according to phase of growth: exponential (FE/DE) or stationary (FS/DS). See Table SIA and B in the supplemental material for the complete set of data. The gene is 
overexpressed in EE or FS if the ratio is above 1. 

Determined by Fisher exact test, in which statistically significant differential expression between two conditions was determined by looking at genes with FDRs of <0.05 (see 
Materials and Methods for details). 

— , gene with a ratio of < 1 or with no significantly increased expression under those conditions. 
" TCS, two-component regulatory system. 



mer deposition, reverse TCA, or nitrate reductase activity, among 
others (35, 37). The induced expression of 11 genes with allegedly 
reductive sink activity, 8 of which are involved in polymer depo- 
sition, might indicate the participation of reductive stress during 
growth of M. tuberculosis in LC-FA (Table 3). Our results suggest 
that polymer deposition, which includes TAG accumulation, as 
well as the formation of other complex lipids, such as SL-1, PAT, 
DAT, and PDIM (Table 3), was the main reductive sink mecha- 
nism used by M. tuberculosis when grown on fatty acids as carbon 
sources. This is consistent with the observed lipid droplets inside 
M. tuberculosis during the stationary phase of growth in LC-FA 
medium (Fig. 3A) and agrees with the suggestion that lipid accu- 
mulation could be a reputed marker of M. tuberculosis during 
latent infection (17, 38). The accumulation of lipid droplets inside 
M. tuberculosis in our LC-FA in vitro model, together with the high 
expression of tgs-1, prompted us to test for drug tolerance, an- 
other dormancy-related characteristic of the bacilli. To do this, we 
determined the drug susceptibility of M. tuberculosis to a mixture 
of four drugs known to be active against replicating and nonrep- 
licating bacilli (39). Our results showed that M. tuberculosis man- 
ifested greater drug tolerance in the presence of lipids than during 



growth in dextrose (Fig. 3B and 3C). The relationship between the 
detected overexpression of tRNAs and drug tolerance deserves 
further investigation, as it represents an attractive idea regarding 
the participation of tRNAs in the dormant stage of M. tuberculosis. 

Comparison of the LC-FA model with other dormancy- 
related in vitro models. The host environmental cues that trigger 
the M. tuberculosis dormant stage or long-term infection are not 
fully elucidated, although changes associated with hypoxia, star- 
vation, and low pH have been shown to be important ( 1 7 ) . To test 
the suitability of growth in LC-FA as an in vitro model, we com- 
pared the gene expression profile of our lipid-enriched condition 
with other well-established in vitro models designed to mimic the 
hypoxic and starvation conditions faced by the pathogen within 
the host. Genes with significantly higher expression during growth 
in LC-FAs (see Table SIE and F in the supplemental material) 
were compared with the highest-scoring genes derived from array 
data obtained from the meta-analysis published by Murphy and 
Brown (40). The total number of genes shared by any two condi- 
tions was very similar, with no significant differences in percent- 
ages (see Fig. S4A). Therefore, growth in LC-FA can be another 
suitable in vitro model for the analysis of genes related to the dor- 
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TABLE 3 Genes with significantly increased expression and their putative relation to reductive stress 





Gene 


Growth 


Functional 




Putative influence 


RPKM 




Gene tag 


name 


phase 


category* 


Function 


on reductive stress"^ 


ratio'' 


P value' 


Rv0243'' 


fadA2 


FS 


1 


Lipid degradation 




1.60 


9.87 E-05 


1\ V / w 




FS 


2 


T3inr»riii(~lpncinp-ninnr\<;nnatp i*pnii(~ta*;p 


NADH/NADPHX 


2.67 


0.004 






FE 


7 


T — T Oi — hci+P H P n \/H VflfTPTl d CP 




1.09 


0.029 


Rv0697 




FE 


7 


Unknown; probable dehydrogenase 




4.33 


0.001 


Rvl 1 RD 

IVV 1 1 ou 


pks-3 


FE 


I 


TJi"^^r\a faIp I \?I/'p+i H p cimtncicp 
I^lUUdUlC UUiyKCLlQC aVllLlldSC 


P/^K/TYIPf H PT^/^ Cl1"l /^n fcMd 
Uiy lllCl LiCUUSlLlUllMP 


3.20 


5.9 E-04 






FE 


7 


Probable GADPH^ 


NADH^I^ 


1.15 


0.028 


Rvl652'' 


(ITgC 


FS 


7 


Argininc biosynthesis 


N AD H/N ADP H(|N 


3.00 


8.33 E-05 


Rvl QflRr 

XVV ± 


katG 


FE 


0 


(^atal acp/npm Yin ^icf 




1.57 


0.004 


Rv2382c 




FE 


1 


Polylcetide synthsse 


Polymer depositions^ 


4.00 


0.032 


Rv2590 


fadD9 


FE 


I 


IjlLJiU UCgi dUd-liUll 


MAHH/NADPHA 


3.71 


5 68 F-n7 


Rv2781 




FE 


7 






3.80 


1.94 E-10 


Rv2794c 


nritT 
ppii 


FE 


I 


4 -Pnncnnr»r>antptnpin vl tranctpfacp 

t: JTlUJoLJll^JL/aJlLCLllClliyi LI alliMCl tijC 


Pnl^/mpf npnncitir»n fcM< 
IT uiyiii^i UC IJU&l LlWll 


1.18 


0.041 






FS 








3.00 


8.33 E-05 


Rv2931 


ppsA 


FS 


I 


PT~)TNA nincvTitnpci<; 

IT ivllVl UlUO y 11 lllColo 


Pnli/mpi" npnncitinn kM< 

IT UiyillCl Lie IJUal LlWli Hp 


2.00 


0.041 


Rv2995c 


leuB 


FS 


7 


Leucine synthesis 


NADH/NADPH(|> 


2.60 


6.54 E-05 


Rv2996c 


serAl 


FS 


7 


L-Serine synthesis 


NADH/NADPH(|» 


3.00 


0.019 


Rv3130c« 


tgs-1 


FS 


1 


TAG synthase 


Polymer deposition«|« 


2.31 


2.07 E-13 


Rv3131« 




FE 


10 


Unknown 


Reductant collector*!* 


2.89 


1.09 E-12 






FS 




Putative NAD(P}H nitroreductase 




4.00 


9.04 E-13 


Rv3145 


nuoA 


FE 


7 


Probable NAD H- dehydrogenase 


Reductant collector(|( 


1.03 


0.008 


Rv3229c 


desA3 


FE 


1 


Lipid desaturase 


Polymer deposition(|( 


3.80 


1.09 E-12 


Rv3230 




FE 


7 


Oxidoreductase 




1.32 


0.006 


Rv3371 




FS 


1 


TAG synthesis 


Polymer deposition«|( 


2.00 


0.041 


Rv3825c 


pks-2 


FE 


1 


Probable polylcetide synthase 


Polymer depositioni|« 


1.27 


1.34 E-04 



" Member of the DosR regulon. 

Functional categories based on http://www.tuberculist.epfl.ch: 0, virulence, detoxification, and adaptation; 1, lipid metabolism; 2, information pathways; 7, intermediary metabo- 
lism and respiration; 10, conserved hypothetical proteins. 

RPKM ratio according to phase of growth: exponential (FE/DE) or stationary (FS/DS). The gene is overexpressed in FE or FS if the ratio is >1. See Table SIA and B in the 
supplemental material for the complete set of data. 

Arrows indicate the putative contributions of genes to the increase (up arrow) or decrease (down arrow) in reductive equivalents. Absence of an arrow indicates a gene that could 
contribute to both an increase and a decrease in reductive equivalents, for example, an oxidoreductase-encoding gene. 

^ Determined by Fisher exact test, in which statistically significant differential expression between two conditions was determined by looking at genes with FDRs of <0.05 (see 

Materials and Methods for details). 

/GAPDH, glyceraldehyde 3-phosphate dehydrogenase. 



mant stage of M. tuberculosis. It is noteworthy that growth in LC- 
FAs and growth under hypoxia conditions shared more genes be- 
longing to the DosR regulon than either condition compared to 
growth during starvation (see Fig. S4A). Interestingly, one of the 
LC-FAs used was palmitate, a fatty acid present in the culture 
medium used in one of the classical in vitro models of persistence 
(41). 

Because of the suggested involvement of reductive stress in the 
growth of M. tuberculosis in a lipid-enriched medium, we sought 
to identify genes related to reductive stress in the same hypoxia 
and starvation data set used for the above analysis (40). Of the 371 
genes associated with reductive stress in the M. tuberculosis ge- 
nome (http://www.tuberculist.epfl.ch), we excluded genes be- 
longing to the functional category of "lipid metabolism" to avoid 
bias toward the fatty acid model. We compared the remaining 216 
genes with the 206 highest-scoring genes for all three data sets (see 
Fig. S4B). The number of genes related to reductive stress identi- 
fied in each case was different, with the highest represented by 
LC-FA (24/206 genes, 11.1%), followed by hypoxia (14/206 genes, 
6.5%), and the lowest represented by starvation (7/206 genes, 
3.2%). The differences among percentages were close to signifi- 
cant between hypoxia and LC-FA and clearly significant between 
starvation and LC-FA (see Fig. S4B). Recently, it has been shown 
that reductive stress could play a role during long-term mycobac- 
terial infection and that lipids could be stakeholders in this process 



(37). Our results also indicate that growth in LC-FAs might be a 
more suitable in vitro model to better test the putative involve- 
ment of reductive stress during the long-term infection of M. tu- 
berculosis. 

Genes involved in the adaptation of M. tuberculosis to a lipid 
environment. The global analysis of the M. tuberculosis transcrip- 
tome showed increased expression of genes involved in various 
metabolic processes and identified specific genes that provide in- 
sight into cellular adaptation during growth with fatty acids as the 
sole carbon sources (Fig. 2; see Table SIE and F). Searching for a 
set of genes related to the pathways used by M. tuberculosis to 
adapt to a lipid environment, we identified 14 genes that had sig- 
nificantly increased expression in FE and FS than in DE and DS, 
respectively (Table 4) (42, 43). These genes therefore character- 
ized the adaptation of the bacilli to growth in LC-FA medium and 
thus provide a lipid signature for M. tuberculosis; this overexpres- 
sion was confirmed by quantitative reverse transcription (qRT)- 
PCR (see Fig. 55). 

Several of these lipid signature genes are involved in cellular 
intermediary metabolism. This includes the genes pckA and id, 
indicating that bacteria are actively using the glyoxylate shunt, 
consistent with the central role of these enzymes in M. tuberculosis 
lipid metabolism (13, 32). The genepckA is required for growth in 
fatty acids, and mutants defective in pckA are attenuated in both 
macrophage and mouse infections (32). Moreover, id and pckA 
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FIG 3 Dormancy-related phenotype of M. tuberculosis grown in the fatty-acid model. (A) Electron micrographs of M. tuberculosis grown under different carbon 
sources. White arrows indicate presumed lipid bodies identified in LC-FA at stationary phase. Magnification: X60,000 (DE, DS, and FE) or X 12,000 (FS). 
Exponential- and stationary-phase cultures were incubated with drugs in the presence of dextrose or fatty acids. Cultures were collected after 7, 14, or 21 days of 
exposure to drugs as estimated by CFU counting (B) or by calculating the MPN/ml (C) (see Materials and Methods for details). The R, MX, MZ, and AK 
concentrations used were 8, 4, 8, and 8 fig/ml, respectively. Means and standard deviations of four experiments are shown. 



are also upregulated in nongrowing bacilli in murine models (13, 
15). 

Two lipid signature genes belong to the functional category of 
"lipid metabolism" (http://www.tuberculist.epfl.ch): Rvll84, 
which codes for an exported protein that is able to perform lipid 
degradation and pptT, which encodes a phosphopantetheinyl 
transferase that activates the polyketide synthases (44). Two addi- 
tional genes are related to membrane activity allegedly involved in 
either the transport of metabolites (Rvl 733) or associated with the 
ESX-1 secretion system of M. tuberculosis (Rv3614c), the latter 
coding for EspD. 

WhiB3 was the single transcriptional regulator detected as part 
of the signature. This heme-based sensor, which is associated with 
intracellular redox balance, links environmental signaling with 
cellular metabolism (34). Interestingly, five of the members of the 
lipid signature belong to the DosR hypoxia regulon. This suggests 
that induction of the DosR transcriptional regulator could result 
from cotranscription from the upstream Rv3134c gene (45, 46), 
which is part of the lipid signature. The Rv3134c gene harbors a 
usp domain, as does another gene of the lipid signature, Rvl 996. 
One of the DosR-regulated genes, Rv0079c, has recently been re- 
ported as a dormancy-associated translation inhibitor (DATIN) 



that inhibits in vitro protein synthesis. Overexpression of a recom- 
binant DATIN in Escherichia coli and in Mycobacterium bovis BCG 
decreased the growth rate of these microorganisms (42). We 
therefore suggest that the overexpression of this gene in both the 
exponential and stationary phases could be the responsible for the 
observed slower growth of M. tuberculosis in LC-FA (see Fig. SI in 
the supplemental material). Many of the genes included in the 
lipid signature are of unknown function, showing our scarce 
knowledge of how M. tuberculosis adapts to a lipid-enriched con- 
dition. One of these genes, Rvl 733c, codes for a hypothetical pro- 
tein reported to be induced by hypoxia and starvation (5). 

To investigate whether these genes could be coregulated at the 
transcriptional level, we searched for a possible consensus motif 
500 nucleotides (nt) upstream of each open reading frame (ORF) 
(47). We identified a DNA motif in 12 of the 14 genes analyzed 
(Fig. 4). The probability that this motif could be present in these 
genes by chance was lE-07, thus suggesting a biological function 
of the motif in genes of the lipid signature. Interestingly, the motif 
was not found upstream of whiB3, opening the possibility that this 
gene product could be involved in transcriptional regulation of 
the lipid signature genes. The identification of an antisense tran- 
script of the whiBS gene in DS (Fig. 4B.2) strongly suggests the 
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TABLE 4 Lipid 


signature of M. tuberculosis i 


grown in the fatty Ecid model 












Gene 


Functional 




RPKM ratio' 




P value"* 












VJCIIC Lug 


name 


category'' 


Pi in r\Ti 1 rf^T^ff^T^ t^f^ \ 
JTUllUlUJll \_1C1C1C11UC/ 


FE/DE 


FS/DS 


PF v; DF 


FS v; n9 


RvOOTg" 




10 


Unknown; probable transduction 


1.59 


2.10 


01.09 E-12 


3.37 E-07 








regulatory protein (42) 










Rv0211 


pckA 


7 


Phosphoenolpyruvate carboxykinase 


1.03 


1.83 


0.021 


4.03 E-14 


Rv0467 


id-1 


7 


Isocitrate lyase 


3.88 


1.52 


1.58 E-08 


3.89 E-09 


Rvl066 




10 


Unknown; rhodanase like domain"^ 


2.11 


4.10 


1.28 E-07 


9.04 E-13 


Rv 1184c 




3 


Unknown; possible exported 


4.80 


5.00 


3.69 E-07 


0.024 








protein, probable lipase, esterase 
















or quitinase (43) 










Rvl733c« 




3 


Unknown; probable conserved 


2.77 


1.50 


2.62 E-08 


0.012 








transmembrane protein 










Rvigge" 




0 


Unknown; universal stress protein 


1.93 


1.19 


1.09 E-12 


6.27 E-06 








family 










Rv2794c 


pptT 


1 


4 ' -Phosphopantetheinyl transferase 


1.18 


3.00 


0.041 


8.33 E-05 


Rv3131« 




10 


Unknown; putative NAD(P)H 


2.89 


4.00 


1.09 E-12 


9.04 E-13 








nitroreductase 










Rv3134c'' 




10 


Unknown; universal stress protein 


1.03 


3.43 


5.99 E-06 


3.19 E-10 








family 










Rv3416 


whiB3 


9 


Redox sensor; transcriptional 


2.50 


1.09 


0.003 


0.009 








regulatory protein 










Rv3612c 




10 


Unknown 


2.00 


1.10 


0.029 


0.014 


Rv3614c 


espD 


3 


Unknown; ESX- 1 secretion 


2.81 


1.50 


1.09 E-12 


9.04 E-13 








associated protein 










Rv3686c 




10 


Unknown 


1.82 


3.00 


1.60 E-08 


3.26 E-07 



" Member of the DosR regulon. 

^ Functional categories based on http://www.tuberculist.epfl.ch: 0, virulence, detoxification, and adaptation; 1, lipid metabolism; 3, cell wall and cell processes; 7, intermediary me- 
tabolism and respiration; 9, regulatory proteins; 10, conserved hypothetical proteins. 

RPKM ratio according to phase of growth: exponential (FE/DE) or stationary (FS/DS). The gene is overexpressed in FE or FS if the ratio is >1. See Table SIA and B in the 
supplemental material for the complete set of data. 

** Determined by Fisher exact test, in which statistically significant differential expression between two conditions was determined by looking at genes with FDRs of <0.05 (see 

Materials and Methods for details}. 

" http://genome.tbdb.org, accessed December 2012. 



involvement of noncoding RNA in whiBS gene regulation under 
lipid-enriched conditions. 

To assess the relevance of these genes for M. tuberculosis 
growth, we next compared this list with previously reported stud- 
ies aimed at identifying essential genes. Interestingly, none of the 
1 4 lipid signature genes were found to contain IS6110 insertions in 
the in vivo survey of 533 isolates worldwide, with the exception of 
Rvl733, which had an IS6110 insertion in 2 of the 533 isolates 
analyzed (48). In contrast, 10 of the 14 genes are considered non- 
essential for in vitro M. tuberculosis growth, based on transposon 
mutagenesis (49). These data support our proposal that the 
LC-FA model could be more closely related to the conditions that 
the bacilli confront in vivo during TB infection. 

Final remarks. On the basis of the results presented here, we 
propose the following scenario for M. tuberculosis growing on 
even-length LC-FAs (Fig. 5). Lipids surrounding the bacilli, TAG 
or LC-FAs, are hydrolyzed by lipases (such as Rvll84) and then 
internalized through receptors (Rvl733c). The lipids are subse- 
quently converted to acetyl-CoA by j3-oxidation and consumed 
by the glyoxylate cycle (IcU). At this point, PckA can regenerate 
phosphoenolpyruvate to continue the cycle. As a consequence, 
both CO2 and the reducing equivalents NADH and NADPH in- 
crease and are sensed by the DosR, RvOOB 1 , and WhiB3 transcrip- 
tional regulators. These and probably other stress conditions in- 
duce the expression of genes such as pptT, pks-2, pks-3, and tgs-1, 
which are responsible for the synthesis of complex lipids that bal- 
ance the redox state of the cell by polymer deposition. The detec- 



tion of espD suggests secretion, possibly of toxic compounds. On 
the other hand, several genes of the DosR regulon are induced: 
stress proteins (Rvl996, and Rv3134c, which are cotranscribed 
with dosR) and other reductive sinks (Rv3131), inhibitors of pro- 
tein synthesis (Rv0079), and heme stabilization proteins (fdxA). 
Proteins required for porphyrin synthesis (HemA and HemC) are 
also induced (see Table S4B). Finally, the high expression of 
tRNAs could modify membrane properties, resulting in tolerance 
to environmental stresses and antibiotics. We propose that when 
M. tuberculosis grows on LC-FAs as sole carbon sources, the cells 
change their gene expression profile by inducing genes required 
for stress protection, as well as genes required for its survival dur- 
ing long-term infection. Our results therefore provide insight into 
the physiological responses of tubercle bacilli when grown in a 
lipid-rich environment, a poorly studied aspect of TB infection. 

Studying the metabolic activity of M. tuberculosis during latent 
infection is not an easy task (50). TB was recently described as a 
continuum, from active disease to latent infection, and the accu- 
mulation of data indicates that latency might not be considered a 
single condition either (50). Although it is still unclear, it appears 
that dormant bacilli include cells with different metabolic activi- 
ties, from slowly growing to nondividing stages that vary in terms 
of redox balance and require responses to oxidative stress or re- 
ductive stress, balanced by the cooperation of main regulators, 
such as WhiB3, DosR, and the recently identified RvOOBl (6). The 
development of new drugs against dormant bacilli has been diffi- 
cult probably because of our poor knowledge of real metabolic 
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FIG 4 Features of lipid signature genes. (A) Sequence logo of tlie lipid signature motif Ttie nucleotide sequence was determined by performing Gibbs sampling 
of the 500 bp upstream of the initiation codon of the lipid signature genes (see Materials and Methods for details). In the representation of the consensus sequence 
shown in the graph, the height of each nucleotide indicates the level of conservation. (B) Differential transcription of the whiB3 gene. Artemis view of the reads 
mapped to the whiB3 gene under the conditions tested. The blue arrow indicates the orientation of whiB3 transcription. (B. 1 ) Exponential phase: red, dextrose; 
purple, LC-FA. (B.2) Stationary phase: green, dextrose; blue, LC-FA. Antisense transcription is identified in the stationary phase of growth in dextrose. 



activity of M. tuberculosis during the complex period of latency, sent a novel tool suitable for the study of the susceptibility of 
The data derived from the use of the LC-FA in vitro model repre- M. tuberculosis to drugs during latency from a new perspective. 




FIG 5 Proposed scenario of M. tuberculosis growth in the fatty acid model. Growth of M. tuberculosis during the in vitro model of even-LC-FA is shown. See the 
text (final remarks) for details. 
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In summary, our results show a new scenario for the adapta- 
tion of M. tuberculosis to a lipid-rich environment, presenting 
evidence of the role of reductive stress and giving insights into the 
participation of novel molecules such as noncoding RNAs and 
tRNAs. These could play important roles in the adaptation of the 
bacteria for long-term survival. During the coevolution of the hu- 
man host and the tubercle bacilli, these adaptations have probably 
been crucial for bacterial survival and the establishment of a bal- 
anced coexistence, represented by a successful infection without 
disease, which we know as latent TB. 

MATERIALS AND METHODS 

Bacterial strain and culture conditions. M. tuberculosis H37Rv (donated 
by The Pasteur Institute, Paris, France) was grown at 37°C in 100 ml of 
Dubos broth (Difco) supplemented with 10% ADC (Middlebrook) in a 
1 -liter flask to an optical density at 600 nm (ODg^o) of 0.5. The culture was 
divided in two, and cells were pelleted by centrifugation and resuspended 
in (i) Dubos broth base (without glycerol) (BD Difco) containing 0.5% 
albumin (Sigma Fraction V, fatty acid free) and supplemented with a 
mixture of 0.003% FA with even an number of carbon atoms (0.001% 
each oleic [18 carbon atoms], palmitic [16 carbon atoms], and stearic [18 
carbon atoms] acids) or (ii) Dubos medium with 0.5% albumin- 0.75% 
dextrose- 0.04% catalase (ADC enrichment). Cultures were grown with 
agitation at 200 rpm at 37°C, and growth was monitored by measuring the 
OD,„o. 

RNA isolation. Total RNA was isolated from cultures in the exponen- 
tial and stationary phases as previously described (51). Basically, cultures 
were harvested by centrifugation and pellets were resuspended in guani- 
dinium chloride buffer (6 M guanidinium chloride, 0.1% Tween 80, 
1 mM 2-mercaptoethanol, 10 mM EDTA) in a proportion of 1 ml of 
buffer/ 100 ml of culture. Cells were lysed mechanically in a FastPrep 
(Thermo Scientific) with 150- to 200-/L(,m glass beads (Sigma-Aldrich) by 
performing four lysis cycles of 15 s each at high speed (6.5 m/s). Nucleic 
acids were purified with phenol-chloroform-isoamyl alcohol (25:24:1), 
and RNA was differentially precipitated with 0.4 volume of absolute eth- 
anol (added drop by drop). Finally, RNA was washed three times with 
Trizol reagent (Invitrogen). RNA integrity was analyzed with a bioana- 
lyzer (Agilent Technologies) and quantified by spectrophotometry with 
the NanoDrop ND-1000 (Thermo Scientific). 

Library preparation and sequencing. Strand-specific RNA-seq li- 
braries were prepared in accordance with the protocol described by Wald- 
bauer and coworkers (52). Briefly, The DNase-treated RNA (100 ng) was 
fragmented by divalent-cation hydrolysis (Fragmentation Buffer; Am- 
bion, Austin, TX) at 70°C for 12 min to yield fragment sizes between 50 
and 300 nt. After precipitation with ethanol, fragments were subjected to 
poly(A) tailing and end repairing (NEB Reagents). RNA was treated with 
Antarctic phosphatase (New England Biolabs) and then phosphorylated 
at the 5' RNA end with T4 polynucleotide kinase. The transcribing strand 
was labeled by ligation of a 5 ' hybrid DNA-RNA primer and after purifi- 
cation with RNA clean XP beads (Beckman Coulter Genomics), the first- 
strand cDNA synthesis reaction was carried out with Superscript II re- 
verse transcriptase (Invitrogen) and lUumina's poly(T) primer and 
deoxynucleoside triphosphates (20 mM). Reaction components were re- 
moved with Agencourt AMPure XP SPRI beads (Beckman), and primary 
transcripts were enriched with Pfu Hi Fidelity polymerase (Invitrogen) 
and the lUumina spacers as primers. lUumina adaptors and bar codes were 
ligated by PCR in accordance with the manufacturer's instructions. Li- 
braries were purified with SPRI beads (Beckman) and quantified with a 
NanoDrop (Bio-Rad). Sequencing was performed with the lUumina tech- 
nology at Ambry Genetics. 

RNA-seq data analysis. Sequence reads were processed to remove 
poly (A) and spacer sequences, and only those with quality scores of >30 
(lUumina 1.3 + ) and a minimal length of 50 bases were used for analysis. 
Assessment of read quality was done with the FASTX toolkit v. 0.0.13 
(http://hannonlab.cshl.edu/fastx_toolkit/index.html). The reads were 



mapped against the M. tuberculosis H37Rv reference genome (accession 
no. NC000962.2) with Bow tie v. 0.12.7 (53) and visualized with the Ar- 
temis sequence visualization and annotation tool (54). Reads that mapped 
to more than one site where excluded. For the analysis, only those reads 
for which 50% of the sequence length fell within the annotated ORE were 
considered part of the ORE. For the identification of novel expressed IGRs 
in the transcriptome, clean reads were analyzed with an in-house algo- 
rithm that identifies genomic regions with a minimum depth of 10 reads 
and a minimum length of 50 nt. These zones were then annotated by 
comparison with the reference genome, with annotation coming from 
literature and the TubercuList database (http://www.tuberculist.epfl.ch). 
For example, transcripts are annotated as novel intergenic if they show 
expression according to the criteria mentioned and they do not overlap 
any previously annotated feature in the genome. 

Measurement of drug activity. The activity of a drug mixture contain- 
ing rifampin (R), moxifloxacin (MX), metronidazole (MZ) and amikacin 
( AK) against exponential- and stationary-phase cultures of M. tuberculosis 
grown in dextrose or FA medium was tested. R, MZ, and AK were pur- 
chased from Sigma Chemicals (United States); MX was obtained from 
Bayer (Mexico). R was dissolved in ethanol; MZ, AK, and MX were dis- 
solved in distilled water. To determine drug activity, 4- to 6-day-old 
(exponential-phase) cultures and 15-day-old (stationary-phase) cultures 
(containing at least 1X10'^ CEU/ml) were incubated with the drug com- 
bination for 0, 7, 14, or 21 days; R, MX, MZ, and AK were used at their 
maximum drug concentration in serum at 8, 4, 8, and 8 /J.g/ml, respec- 
tively (39). After incubation, culture CFU counts per milliliter were de- 
termined by plating appropriate dilutions as described previously (55). In 
parallel, 1 ml of each culture was taken to perform 10^' to 10^^ dilutions 
with Dubos broth. One miUiliter of every dUution was transferred to each 
of three tubes containing 9 ml of Dubos broth with or without the fatty 
acid mixture, accordingly. These broth cultures were examined for visible 
turbidity every 7 days over a 3-week period of incubation at 37''C. The 
most probable number (MPN) of viable bacilli was then estimated while 
taking in account the patterns of positive and negative tubes (http://www 
.fda.gov/food/foodscienceresearch/laboratorymethods/ucml09656.htm; 
accessed lanuary 2013). Plating of different culture samples on blood agar 
medium showed the absence of bacteria other than mycobacteria. Growth 
of M. tuberculosis in turbid tubes was confirmed by Ziehl-Neelsen staining 
and colonial morphology on Middlebrook 7H1 1 

Analysis of the upstream regions of the genes belonged to the lipid 
signature. For each gene, 500 nt were taken upstream of the initiation 
codon. These regions were fed into the Gibbs sampling algorithm, and the 
sequences were analyzed to detect any possible bias in their composition. 
This was done by means of the program "unifiedcpp." Once corrected for 
possible bias, the samples were run with the following parameters: k-mer, 
12 to 24 nt; -n, for nucleic acid alphabet; -E 3, maximum number of sites 
to find; -S 20, number of seeds to try; -i 1,000, number of iterations. The 
probability of the motif arising by chance was determined by means of the 
hypergeometric distribution (56). For this, the number of genes in which 
the motif was present in the 500-nt upstream region in H37Rv was calcu- 
lated. It turned out that 767 of the 4,060 genes had the motif (this is the 
population distribution), so we wanted to calculate the probability of 
drawing 12 out of 14 genes given the population distribution (767/4,060). 

qRT-PCR of selected genes. The total transcripts of selected genes 
were measured by real-time qRT-PCR. Quantification was performed 
with gene-specific primers (see Table S5) and SYBR green (Molecular 
Probes, Inc.). Samples were subjected to 40 cycles of amplification (dena- 
turation at 95''C for 30 s, specific annealing temperature for 15 s, and 
extension at 72''C for 15 s) with a final extension at 72°C for 5 min. To 
ensure that the fluorescence levels detected were due to the amplification 
of a specific product, a melting curve analysis was performed. Absolute 
analysis was carried out. To do this, a standard curve was obtained for each 
set of primers by using 10-fold dilutions of known amounts of M. tuber- 
ct//osis H37RVchromosomalDNA(l,000, 10,000, 100,000, and 1,000,000 
theoretical copies). Threshold cycle values of each RNA quantitation were 
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interpolated to standard curve to obtain gene expression (number of gene 
copies per microliter). Normalization of these data was performed by 
using 16S rRNA expression levels. 

Transmission electron microscopy. Cells were harvested by centrif- 
ugation at 16,000 X g, washed in sterile phosphate-buffered saline (PBS), 
and fixed in 2.5% (vol/vol) glutaraldehyde- 0. 1 M sodium phosphate for 
1 h. Cells were recovered and resuspended in 1% (wt/vol) osmium tetrox- 
ide for 1 h. Cells were then dehydrated through treatment with a series of 
ethanol solutions (30, 40, 50, 60, 70, 80, and 90%); each treatment was 
performed twice for 10 min each time. A final dehydration step was done 
twice with 100% ethanol for 15 min. Samples were embedded in EPON 
812 resin (EMS). This resin was polymerized at 60°C for 24 h. Resin blocks 
were cut in a Leica Ultracut UCT Ultratome to a 70-nm thickness, placed 
on 200-mesh copper grids (EMS), and stained with 30% (wt/vol) uranyl 
acetate-70% (vol/vol) methanol and counterstained with Reynolds lead 
citrate buffer. Samples were examined with a JEOL JEM- 1010 transmis- 
sion electron microscope with an acceleration voltage of 60 kV. 

Statistical analysis. Statistical analysis for differential gene expression 
was evaluated by using RPKMs as expression values and Fisher's exact test 
on a two-by-two contingency table (57) on a per-gene basis in which the 
two columns of the table were separated by reads that overlap the gene's 
ORE and reads that do not overlap the ORE, and the two rows were 
separated by the two compared samples (EE versus DE or ES versus DS). 
Significant differences were determined with Eisher's exact test with the 
additional criteria of a false-discovery rate (FDR) of <5%. The results 
were consistent with those obtained with the DEGseq package for biolog- 
ical replicates (58). 

Nucleotide sequence accession number. The RNA-seq data sets have 
been deposited in the NCBI Gene Expression Omnibus (59) under acces- 
sion number GSE47863. 
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